clear
clc

%% ************** Parameters **************
shape               = 1;
% shape               = 2;
% shape               = 3;
% shape               = 4;

sigma               = 0.001;
eps_r               = 5;
tx1                 = 0.5;
tx2                 = 1.0;
f                   = 2.4e9;

dx                  = 0.05;
ds                  = 0.25;

distance            = 500;
h_max               = 4;
w_max               = 5;
cs                  = func_cs_mapping(shape, dx, h_max, w_max);

k_w0                = 0.9;
bc                  = 1;
tx{1}               = tx1;
tx{2}               = tx2;
tx{3}               = 1;
tx{4}               = 20;

% Rx setup only used for plotting
rx{1}               = tx1;
rx{2}               = tx2;
rx{3}               = 1;
rx{4}               = ds;

%% ************** Simulation **************
[~, rss] = Func_VPE(f, eps_r, sigma, distance, tx, k_w0, shape, bc, cs, dx, ds);

%% ************** Plotting **************
figure;
p_dist = 0:ds:distance;
pl = Func_Extract_Field( ...
    f, ...
    distance, ...
    tx, ...
    rx, ...
    shape, ...
    cs, ...
    dx, ...
    ds, ...
    rss ...
    );

plot(p_dist, pl, '-k', 'LineWidth', 1.5);

grid on
xlim([0 distance]); ylim([-80 0]);
xlabel('Axial Distance [m]')
ylabel('Received Power [dBm]')
pbaspect([1.5 1 1]);

%% ************** Functions **************
function cs = func_cs_mapping(shape, dx, h_max, w_max)
    if shape        == 1
        height      = h_max;            % height of the tunnel (m)
        width       = w_max;            % width of the tunnel (m)    
        cs          = {width, height};

    elseif shape == 4
        x_b1        = 2;                % horizontal distance from first breaking point to the center
        y_b2        = 1.5;              % vertical distance from second breaking point to the ground
        w_half      = w_max/2;          % half-width of the tunnel; if top is semi-circle, set w_half = h_top
        r2_bottom   = w_half;           % radius for the bottom if circular or elliptic 
        h_bottom    = y_b2 - 2*dx;      % height from center to the bottom (platform); if circular, set h_bottom = h_top
        h_top       = h_max - h_bottom; % height from center to the top boundary; if top is semi-circle, h_top is the radius
        cs          = {w_half, h_top, r2_bottom, h_bottom, x_b1, y_b2};
    else
        w_half      = w_max/2;          % half-width of the tunnel; if top is semi-circle, set w_half = h_top
        h_top       = w_half;           % height from center to the top boundary; if top is semi-circle, h_top is the radius
        h_bottom    = h_max - h_top;    % height from center to the bottom (platform); if circular, set h_bottom = h_top
        r2_bottom   = w_half;           % radius for the bottom if circular or elliptic
        cs          = {w_half, h_top, r2_bottom, h_bottom};
    end
end
